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We consider a Kuramoto model for the dynamics of an excitable system consisting of two coupled 
active rotators. Depending on both the coupling strength and the noise, the two rotators can be 
in a synchronized or desynchronized state. The synchronized state of the system is most stable for 
intermediate noise intensity in the sense that the coupling strength required to desynchronize the 
system is maximal at this noise level. We evaluate the phase boundary between synchronized and 
desynchronized states through numerical and analytical calculations. 
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I. INTRODUCTION 

Networks of coupled nonlinear oscillators provide use- 
ful model systems for the study of a variety of phenomena 
in physics and biology [H. Among many others, exam- 
ples from physics include solid-state lasers [2] and coupled 
Josephson junctions |3j,y]. In biology, the central nervous 
system can be described as a complex network of oscil- 
lators [1] , and cultured networks of heart cells are exam- 
ples of biological structures with strong nearest-neighbor 
coupling [6i]. In particular, the emergence of synchrony 
in such networks 0, [1| has received increased attention 
in recent years. 

Disorder and noise in physical systems usually tend 
to destroy spatial and temporal regularity. However, in 
nonlinear systems, often the opposite effect is found and 
intrinsically noisy processes, such as thermal fluctuations 
or mechanically randomized scattering, lead to surpris- 
ingly ordered patterns j^]- For instance, arrays of cou- 
pled oscillators can be synchronized by randomizing the 
phases of their driving forces [l^, Synchronization 
in these systems is caused by the interactions between 
the elements and results in the emergence of collective 
modes. It has been shown to be a fundamental mecha- 
nism of self-organization and structure formation in sys- 
tems of coupled oscillators [l^- Biological systems of 
neurons are subject to different sources of noise, such 
as synaptic noise or channel noise [13] ■ In particu- 
lar, sensory neurons are notoriously noisy. Therefore, the 
question arises how stochastic influences affect the func- 
tioning of biological systems. Especially interesting are 
scenarios in which noise enhances performance. In the 
case of stochastic resonance fl3], e.g., noise can improve 
the ability of a system to transfer information reliably, 
and the presence of this phenomenon in neural systems 
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has been investigated |16l . [17|. Furthermore, numerous 
studies have addressed the effect of noise on the dynam- 
ics of limit cycle systems [H [H [H [H [HI, [H, [I! . 

Small neural circuits composed of two or three neurons 
form the basic feedback mechanisms involved in the regu- 
lation of neural activity [l^] • They can display oscillatory 
activity [1^ [2^ and serve as central pattern generators 
involved in motor control [27j . Here, we consider a sys- 
tem of two limit cycle oscillators with repulsive coupling. 
We investigate the influence of the noise and the coupling 
strength on the dynamics of the system. We distinguish 
between two different classes of dynamics, a synchronized 
state, in which the joint probability density of the oscilla- 
tor phases is characterized by a single-hump shape, and 
a desynchronized state. The single-hump shaped distri- 
bution of the oscillator phases has been modeled by a 
Gaussian distribution [Ij, [ll] , and systems consisting of 
a large number of oscillators were analyzed by examin- 
ing the resulting dynamics for the mean of the oscillator 
phases [2^. In contrast, the simplicity of our two oscilla- 
tor system allows us to obtain the stationary probability 
density function for the full system both numerically and 
analytically. We show that the probability distribution 
of the oscillator phases has the single-hump shape only 
for weak coupling, whereas it deviates from this shape 
for strong coupling. We evaluate the coupling strength 
at which the transition between the two forms of the 
probability distribution occurs as a function of the noise 
intensity. 

In Sect, ini we introduce the Kuramoto model for ex- 
citable systems. Under the influence of noise, the dy- 
namics of the limit cycle oscillators are described by a 
stochastic differential equation (SDK), and we state the 
Fokker-Planck equation for the system. In Sect. IIIIl we 
consider a single active rotator driven by noise and derive 
its mean angular frequency from the stationary solution 
to the Fokker-Planck equation. We compare our ana- 
lytical results with Monte-Carlo simulations of the cor- 
responding SDE. In Sect. IIVI we consider two coupled 
deterministic rotators and perform a bifurcation analysis 
of the system. We show that the system possesses a fixed 



2 



point that is stable for small coupling strengths but looses 
its stability when the coupling is increased. For some 
range of the coupling strength, the stable fixed point and 
a stable limit cycle coexist. In Sect. |Vl we consider two 
coupled active rotators under uncorrelated stochastic in- 
fluences. In Sect. W A\ we solve the Fokker-Planck equa- 
tion of the system numerically and show that the shape 
of the probability distribution undergoes a characteristic 
change, corresponding to the transition from a synchro- 
nized to a desynchronized state, as coupling is increased. 
We evaluate the boundary between the synchronous and 
the asynchronous regime through a Fourier expansion ap- 
proach in Sect. IVBl A summary concludes the paper in 
Sect. IVll 



EXCITABLE SYSTEMS AND THE 
KURAMOTO MODEL 



II. 



Neurons can display a wide range of behavior to differ- 
ent stimuli and numerous models exist to describe neu- 
ronal dynamics. A common feature of both biological 
and model neurons is that sufficiently strong input causes 
them to fire periodically; the neuron displays oscillatory 
activity. For subthreshold inputs, on the other hand, the 
neuron is quiescent. When a subthreshold input is com- 
bined with a noisy input, however, the neuron will be 
pushed above threshold from time to time and fire spikes 
in a stochastic manner. In this regime, the neuron acts 
as an excitable element. In general, an excitable system 
possesses a stable equilibrium point from which it can 
temporarily depart by a large excursion through its phase 
space when it receives a stimulus of sufficient strength 
[23 |. Besides neurons, chemical reactions, lasers, models 
of blood clottingj^ and cardiac tissues all display excitable 
dynamics [H, Ho, HH, [H, . Pulse propagation, spiral 
waves, spatial and temporal chaos, and synchronization 
have been studied in these systems [13, IS ISE [sS] ■ 

The phase dynamics of an active rotator without inter- 
action and random forces can be described by the model 
developed by Kuramoto and coworkers [H, [s^ : 



(f)(t) = uj — asin</)(i) . 



(1) 



To obtain the case of the excitable system with one 
stationary point, one chooses the parameter a > u. 
When we have n coupled identical oscillators, subject 
to stochastic influences, the model is described by the 
Langevin equation [23| 



4>i{t) = w — asin0i(i) 



(2) 



Here, we take the rji to be uncorrelated Gaussian white 
noise, i.e., {tuit)) = 0, (r?i(ti)77j (^2)) = 2cr^(ti - t2)5ij. 
We will concentrate on the simplest case, namely that the 
coupling functions Wij are sin-functions multiplied by a 
coupling constant Wij, i.e., Wij(0) = WijSa\(j}. Then, 



the dynamical evolution of the system's probability den- 
sity function P{4>, t) is described by the Fokker-Planck 
equation 

= (3) 

i=i 

n n „2 

+EEa^[A.(</')mt)], 

where in our case the drift terms read 

n 

Di{(j)) = UJ - a sin cj)., - ^ sin[cj)j - cj).i) (4) 



and the diffusion terms are given by 



(5) 



Since the angle variables describe the phases of the 
oscillators, the probability density function must satisfy 
the periodic boundary conditions 

P(0i, ••• , 0, = 0, ••• , (6) 
= ^"(01, ■ ■ • , '/'i = 27r, • • • , 0„, i) , i = 1, • • • , n . 

Furthermore, the normalization condition for the proba- 
bility density reads 



Jo Jo 



III. SINGLE-ROTATOR SYSTEM 



(7) 



We first exam a single rotator subject to a noisy input 
and, following Ref. [40|, calculate the mean frequency of 
oscillations as a function of the noise level. In this case, 
the Fokker-Planck equation ^ reads 



d_ 
di 

with 



1 

[Di(b)P{(b,t)]+a—Pi<j,,t), (8) 



— Lo ~ a sin ( 



(9) 



We can thus write the drift term as the negative gradient 
of a potential, D = —dV/dcj), with the potential given by 

V{(f)) ^ ~uo(t)~ a cos (j) + c. (10) 

Introducing the probability current 

5(0,t) = i?(0)P(0,t)-a^P(0,i), (11) 

the Fokker-Planck equation takes the form of a continuity 
equation. 



|mt)-fA5(^,t) = o. 



(12) 
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We now look for a stationary solution of the form 
P((^, t) = P{4>), S{(j), t) — S{(f>). In this case, we conclude 
from (fT2| that the derivative of the probability current 
with respect to (f) must vanish, and we have to solve 



(13) 



The constant probability current S is related to the mean 
drift velocity, i.e., the mean angular frequency of the ac- 
tive rotator system according to w = 2ttS. The solution 
to the ordinary differential equation (|13p is given by 



P(0) = Ce- 



S 



e ' . (14) 



The integration constant in (|TOl) can thus be absorbed 
into the constant C in (fT4|) . and the two free constants 
S and C are determined by the periodicity and normal- 
ization conditions ^ and ([7]). These two conditions can 
be written in matrix form as 



j„ d0e - Jo d0Jo d(j)e ^ H e 
: " — e " acp e 



.(15) 

Denoting the determinant of the 2x2 matrix in the last 
expression as det, the constants C and 5* are given by 



C 

s 



det 

(T 

dot 



27T 



(16) 
(17) 



Speciahzing to the potential of the active rotator (|T0)) . 
we obtain 



27rcr 



27r 



s((?i)+^') — cos 0] 



(18) 



Note that in the limit ct — + oo the integrand in the de- 
nominator approaches one, and Co converges to cj. To ob- 
tain the leading order behavior of u) in the limit of small 
noise, we approximate the denominator using Laplace's 
method described in Ref. ^l|. According to Laplace's 
method the asymptotic behavior of the integral 



I{x) = I dtf{t)e''s^'^ 
as a; — > cxD is given by 
/(x)~ 



/2^/(c)e^ff(^) 



(19) 



(20) 



Here, it is assumed that g{t) has a maximum at t = 
c with a < c < b and that /(c) ^ and g"{c) < 0. 
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FIG. 1: Average angular frequency of the single-rotator as a 
function of the noise intensity. The solid line shows the result 
(|18p . The dots represent results from Monte-Carlo simula- 
tions (mean ± standard error of the mean) of the Langevin 
equation ([2| . For each value of the noise intensity, forty runs 
where simulated up to T = 400. The first inset shows a com- 
parison between the asymptotic expansion (|26l dashed line) 
and numerical evaluations of the expression (|18l solid lines) 
for small noise. The diamonds in the second inset show the 
logarithm of the relative deviation between the result (|18p and 
its asymptotic approximation (|26|l . Parameters are: lu — 1, 
a = 1.2. 



Wc first apply Laplace's method to the inner integral in 
the denominator of (US]), which wc denote as I (a). The 
function a[cos{(f) + cj)') ~ cos 0] has a maximum inside the 
interval < g!i < 27r at 



TT + arctan ■ 



sm( 



1 — cos (j)' 
Using ((20|l we thus obtain for cr ^ 

gf [cos(0o+0')-cos(>o]-^<^' 



I{a) 



2-na 



^COs(0o + </>') - COS(/)o 



(21) 



(22) 



The argument of the exponential function in the last 
identity can be simplified to 



a — cos ( 



(23) 



whose maximum within the interval < ^' < Stt is at 

(24) 



2 arccos — 
a 



Using this and applying (pOI to the intermediate result 
(|22)) . we obtain 



/(a) 



\/ a? — Up' 



0. (25) 
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FIG. 2: Stable and unstable fixed points and oscillations in the deterministic two-rotator system, (a) shows the bifurcation 
diagram with stable (solid lines) and unstable (dashed lines) fixed points of the system (|27|) for the choice of parameters a; = 1, 
a = 1.2, Wi2 = W21 ~ w. Dots indicate the minimum and maximum values of oscillations in the value of A that result for 
the initial conditions <1> = 0, A = tt/2. (b) depicts for w = 0.308 the boundaries between the regions in the space of initial 
conditions for which the system converges to the limit cycle or the stable fixed point. 



The leading asymptotic behavior of w as cr 
given by 



^asy 



,2„-|(v/^?^- 



uj arccos - 



is then 



(26) 



Figure [T] shows the mean angular frequency as a func- 
tion of the noise level a. The evaluation of the analytical 
expression (jlSp yields results that are in good agreement 
with Monte-Carlo simulations of the Langevin equation 
(|2|). Furthermore, the asymptotic expansion (|26p is in 
excellent agreement with numerical evaluations of ^ 
for small noise. 



IV. DETERMINISTIC TWO-ROTATOR 
SYSTEM 

We next turn to a system of two coupled active rota- 
tors, where we first consider the deterministic case, i.e, 
CT = 0. In particular, we are interested in rotators with re- 
pulsive coupling, i.e., we consider the case Wi2, W21 > 0. 
Introducing the center of mass and difference coordinates 
$ = (0i + 02)/2 and A = (0i — 02)/2, the set of equations 
([2]) takes the form 

$(t) = Lu - a sin ^{t) cos A{t) 

+ {wi2 - W21) sin A(t) cos A(i) , 

A{t) = -acos^{t)smA{t) 

+ {wi2+W2i)smA{t)cosA{t) . (27) 

The system has a trivial stationary point at <I>(i) = $0 = 
sin~"'^(aj/a), A{t) = 0, whose stability we analyze by lin- 
earizing the system ([27]). Writing ^(t) = $0 + e$(i), 
A{t) = eA{t) we obtain to first order 



dt V eA(t) 

—\J a? — LtP- 




(28) 



W\2 



W12 
' W2I 



W2I 



The real parts of the eigenvalues of the 2x2 matrix 
on the right-hand side of the last identity determine the 
stability of the fixed point ($01 0). Under the assumption 
a > UJ the first eigenvalue Ai — — V — uj"^ is always 
real an d negativ e. The second eigenvalue A2 — W12 + 
W21 — \JaF—uP' is also always real; for small coupling it 
is negative, but when the sum of the coupling strengths 
W\2+W2\ increases it becomes positive and the fixed point 
($0,0) loses its stability in, as it turns out, a subcritical 
pitchfork bifurcation. Further fixed points of the system 
can be determined and turn out to be unstable for all 
values of the coupling strengths. In the case w\2 = W2\ = 
w they are given by 



1 . _i / 'iLJw\ 



Ai 



asin$i 



(29) 



Figure ^a) shows a bifurcation diagram of the system. 
For small coupling strength, the system does not dis- 
play oscillatory behavior. When the coupling strength 
is increased above a critical value, a stable limit cycle 
emerges from a homoclinic orbit. For a small range of 
coupling strengths, the stable fixed point coexists with 
the stable limit cycle. In this case, it depends on the ini- 
tial conditions whether the system will converge toward 
the fixed point ($o,0) or the limit cycle. Figure [2{b) 
shows the attractors for fixed point and limit cycle dy- 
namics in the (<&, A)-plane for W12 — W21 = 0.308. In the 
strong-coupling limit, the minimum and maximum of A 
in Fig.[2Ia) both converge toward 7r/2. Thus, the system 



5 



(a) (b) 




FIG. 3: Synchronized and desynchronized modes in the stochastic two-rotator system. The stationary solution to the Fokker- 
Planck equation @ is shown for different values of the coupling strength. In both (a) and (b), we have ■W12 = W21 = w and 
uj = 1, a — 1.2, a — 0.4. In (a) the coupling strength is w — 0.3 and the rotators are in a synchronized state; in (b) the 
coupling is increased to it; = 0.4 and the two rotators desynchronize. 



approaches antisynchronous oscillatory dynamics, where 
and (f)2 are phase shifted by tt while their sum increases 
constantly. 

V. STOCHASTIC TWO-ROTATOR SYSTEM 

We now consider the coupled two-rotator system in the 
case where both rotators receive uncorrelated stochastic 
driving. The temporal evolution of the probability den- 
sity of this system is given by the Fokker-Planck equation 
([3]) with the drift and diffusion coefficients ([4]) and ([5]). 

A. Numerical Results 

First, we investigate the stationary solution to the 
Fokker-Planck equation numerically. To this end, we nu- 
merically solve the partial differential equation ^ un- 
der the periodic boundary conditions ^ for the homo- 
geneous initial condition P((/)i, (/)2, t = 0) = l/47r^ and 
observe that the solution converges to the stationary so- 
lution after some time. Figure [3] shows the stationary 
solution in the coordinates <& and A for two different val- 
ues of the coupHng strength. We find that, depending 
on the strength of the noise and coupling, two different 
characteristic forms of the stationary solution exist. In 
the case shown in Fig. the probability density is 
peaked around the stable fixed point of the deterministic 
two-rotator system ($o,0). In Fig.[3]Jb), the peak at the 
fixed point ($0j 0) is much less pronounced. Furthermore, 
if we consider the probability distribution for A = ±tt/2, 
i.e., at the edge of the region shown in Fig.[3l we see that 
the probability distribution is not given by one central 
hump anymore. In order to distinguish between the two 
different scenarios in a quantitative way, we consider the 



marginal stationary probability density 

P(A) = / d$P($, A) . (30) 

Figure d] shows this quantity for one level of the noise 
intensity a and for different coupling strengths. For 
weak coupling, P(A) has a pronounced maximum at 
A = 0. For increasing coupling strengths, this maxi- 
mum decreases and eventually turns into a minimum. 
We can thus classify the system dynamics as synchro- 
nized or desynchronized according to the sign of the sec- 
ond derivative of P(A) at the origin and can label the 
a-w plane accordingly. In the next section, we calculate 
the phase boundary between the synchronized and desyn- 
chronized regime through a Fourier expansion approach. 



B. Fourier Expansion Results 

The probability density P(^i, (1)2) is periodic in (jji and 
02, so we expand it as 

P(0i, (/.2) = C(fci, fc2)e^^'=^^^+'=^*^) . (31) 

k-i .k2 

Inserting this approach into the right-hand side of ([3]) 
yields together with ^ and ([5]) 

0= ^ C(/ci,fe)e'(''^'^^+''''^"' (32) 

X |a(cOS0i + C0S(/)2) - {Wi2 + W21) C0S((?!)2 - 4>l) 

— iki[LU — asm(f)i — W12 sin((/)2 — 0i)] 

— ik2[ui — asin(j)2 — W21 sin((/)i — ^2)] — crkf — ak^ > ■ 
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FIG. 4: Marginal probability density for different values of 
the coupling strength w = wi2 — W2i. The coupling strength 
for the curve with the highest value at A = (solid line) is 
It) = 0.1 and increases from curve to curve in increments of 
Sw — 0.2 to the maximum value w = 1.1 (dotted line). Other 
parameters are: uj — 1, a = 0.4, a = 1.2. 



The term inside the curly brackets on the right-hand-side 
of the last identity is itself periodic in 0i and 02 and can 
also be expanded as a Fourier series 

{•••}= E C{h,h,h,k2)e^^''*'+'^'^^^. (33) 

\h\<l.\l2\<l 

Here, the coefficients C{li,l2) read 



2(l±fci), 
|(l±fc2), 



C{0,ki,0,k2) = -iuj{ki + k2) - a{kl + akj) , 
(7(±l,fci,0,fc2) 

C(0,fci,±l,fe) 
C(l,fci,-1,A:2) 
C(-l,A;i,l,A;2) 

C(±l,fci,±l,fc2) 



1 + fci 
2 

1-ki 



Wl2 
Wl2 



l_k2 
2 

1 + A:2 



W21 
W21 



0. 



We can then rewrite ([5^ as 



0- E 



,i(fci0i+fe202) 



(34) 



(35) 



All ,k2 

X ^ C{ki- li,k2- l2)C{li,ki- Ii,l2,k2- h) ■ 

\h\<l,\l2\<l 

Setting the inner sum to zero, we obtain an infinite sys- 
tem of algebraic equations. In order to obtain the A^th 
Fourier order approximation we truncate the outer sum 
such that we set Cjv(fci, ^2) = for \ki\ > N or 1^2] > N. 
Then, we have to solve a system of (27V -|- 1)^ — 1 algebraic 
equations in order to obtain the expansion coefficients to 



A^th order CAr(fci,fc2), where the additional index N in- 
dicates the approximation order. Finally, the coefficient 
Cat (0,0) is determined from the normalization condition 
as CAr(0,0) = l/47r2. 

As an illustrative example we now consider the first 
order in the Fourier expansion for the case 'W12 — W21 = 
w. The system of algebraic equations we need to solve 
then reads 

a[Ci(0,-l) + Ci(-l,0)]-4(<T-icj)Ci(-l,-l) = 0, 
47r2[2(CT-iw)Ci(-l,0) + wCi(0, -1)] = a. 



27r2{a[Ci(-l, 0) + Ci(0, 1)] - 4aCi(-l, 1)} 
47r2[wCi(-l,0) + 2(f7 - iw)Ci(0, -1)] 
47r2[2(f7 + iw)Ci(0,l) + u;Ci(l,0)] 
2^2{a[Ci(G, -1) + C7i(l,0)] - 4aCi(l, -1)]} 



a, 
a. 



47r2[2(cr + icj)Ci(l,0) + u;Ci(0,l)] a, 
a[Ci(0,l) + Ci(l,0)] -4(o• + ^t^)Cl(l,l) = 0. 

(36) 

From this we obtain the first order approximation 



a{2a(3-fa{cos(f>i -I- cos (^2) 



1 

4^ 

-I- 4a/3(Ta;(sin -I- sin(/)2) 
-I- (a^7(T^ - 2a^crw^) cos(0i + 02) 
+ - ■u;/3(7^ -I- 4tj^)] cos((/)i - 

-|- a^a{wuj + Aau)) sin{(l)i + (f>2)} , 



(37) 



with the abbreviations 
1 

" " 47r2a/3(72 +4tj2) ' 



J — w + 2a . 

(38) 



Substituting the coordinates 0i and 02 according to 0i = 
$ -I- A and 02 = $ — A and integrating with respect to 
$ we obtain the marginal probability density 

A (A) ^ ^ + 27ra/3(a^7 - wj^ - Awoj^) cos(2A) . (39) 

Setting its second derivative to zero, we obtain the equa- 
tion 



a^iw + 2a) - w[{w + 2af + 



0, 



(40) 



which we can solve in w or in a. Eventually, we want to 
obtain w as a function of a. However, since we have a 
cubic equation in w and only a quadratic equation in cr, 
for convenience we express c as a function of u: 



- 2w^ ± 



16w2(x)2 



Aw 



(41) 



This procedure can easily be generalized to higher or- 
ders. Figure [5] shows the resulting phase diagram ob- 
tained from solving the Fokker-Planck equation numer- 
ically and from the Fourier expansion. The accuracy of 
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FIG. 5: Regimes of synchronized and desynchronized dy- 
namics. The phase boundary between the synchronized and 
desynchronized regimes is shown as a function of the noise 
strength a. Areas below the curves correspond to the synchro- 
nized, areas above the curve to the desynchronized regime. 
The sohd lines show the results of the first four Fourier orders, 
the dots represent numerical results. The diamond represents 
the coupling strengths for which the fixed point ($o, 0) of the 
deterministic system becomes unstable; the square indicates 
the value of w at which the stable limit cycle is first observed. 
The insets show results for small and for large noise. In the 
first inset (small noise) the results from the fourth and tenth 
Fourier orders are shown. The second inset (large noise) shows 
the results from the first (solid line) and second (dashed line) 
Fourier orders. Parameters are u) = 1, a = 1.2. 



the Fourier expansion results improves with increasing 
strength of the noise. This can be seen, for instance, in 
the second inset of Fig. [51 where even the first expansion 
order yields very accurate results for strong noise. In 
general, even relatively low orders in the expansion give 
a good estimate for the phase boundary for a wide range 
of noise strengths, as can be seen from the results for the 
fourth expansion order in Fig.[5l However, for very small 



noise levels the Fourier expansion diverges, as is exem- 
plified in the first inset in Fig. [5] for the fourth and tenth 
expansion orders. Considering the first inset in Fig. [5l 
we conclude that in the limit ct — > the results from 
the Fourier expansion approach a value of the coupling 
strength for which the stable fixed point coexists with the 
limit cycle in the deterministic system. Therefore, nei- 
ther the existence of the stable limit cycle nor the stabil- 
ity of the fixed point can be used exclusively to determine 
the zero-noise limit of the phase transition between the 
synchronized and desynchronized states. Strong noise 
has a desynchronizing effect on the system, as the mini- 
mal coupling for desynchronization vanishes in the limit 
of cr oo. If the noise is weak, however, it stabilizes 
the synchronized state, as is indicated by the initially 
upward slope of the phase boundary in Fig. [S] In conclu- 
sion, the synchronized state of the system is most stable 
for intermediate noise. 



VI. SUMMARY 

We have investigated the transition from synchronized 
to desynchronized behavior in a system of two-coupled 
active rotators under stochastic influences. The two 
regimes are distinguished by the sign of the second deriva- 
tive of the marginal probability density at vanishing 
phase difference. We have evaluated the phase bound- 
ary between the two states in the (coupling strength) - 
(noise intensity) plane. Finally, we have shown that the 
synchronized state is most stable, in the sense that the 
coupling strength required to desynchronize the system 
is maximal for nonvanishing noise intensity. 
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